#include "kernel/cell.hpp"
#include "yocto/lua/lua-config.hpp"

Cell:: ~Cell() throw()
{
}


#define __luaRead(NAME) Lua::Config::Get<lua_Number>(L,NAME)

Cell:: Cell( lua_State *L ) :
M(0),
nvar(0),
library(),
chemsys(library),
init_in(library),
init_out(library),
solutions(NumSol,as_capacity),
variables(),
enzymes(8,as_capacity),
pAE( new AE(L) ),
EnzAE( pAE ),
pNHE( new NHE(L) ),
EnzNHE( pNHE ),
pNaK( new NaK(L) ),
EnzNaK( pNaK ),
pLactate( new Lactate(L) ),
EnzLactate( pLactate ),

Surface( __luaRead("Surface") ),
Volume(  __luaRead("Volume" ) ),

Em( __luaRead("Em") ),
RT( R * (273.15+25) ),
F_over_RT( Faraday / RT ),
Capa( Surface * __luaRead("Cm") ),
drvs( this, & Cell::computeDiff ),
dt_ctl(0),
dt_max( __luaRead("dt_max") )
{
	
	if( dt_max <= 0 )
		throw exception("Invalid dt_max!");
	
	//-- loading all species
	library.load(L,"library");
	library.close();
	(size_t&)M    = library.size();
	(size_t&)nvar = M + NumExtraVar; 
	variables.make(nvar,0);
	
	//-- loading reactions
	chemsys.load(L,"chemsys");
	chemsys.compile();
	
	//-- loading initializer for inside sol
	init_in.load(L,"init_in");
	init_in.compile();
	
	//-- loading initializer for outside sol
	init_out.load(L,"init_out");
	init_out.compile();
	
	//-- preparing solutions
	for( size_t i=1; i <= NumSol; ++i )
	{
		const Solution sol( library );
		solutions.push_back(sol);
	}
	
	init_in.setup(chemsys, 0.0);
	solutions[ iS_in ].load_from(chemsys);
	std::cerr << "S_in=" << std::endl; solutions[ iS_in ].display( std::cerr );
	
	init_out.setup(chemsys, 0.0);
	solutions[ iS_out ].load_from(chemsys);
	std::cerr << "S_out=" << std::endl; solutions[ iS_out ].display( std::cerr );
	
	//-- build database
	record( EnzAE );
	record( EnzNHE );
	record( EnzNaK );
	record( EnzLactate );
	
	std::cerr << "-- cell Surface=" << Surface << ", Volume=" << Volume << std::endl;
	//exit(1);
		
}

void Cell:: record( const EnzymePtr &enz )
{
	std::cerr << "--> using enzyme '" << enz->name << "'" << std::endl;
	if( !enzymes.insert(enz) )
		throw exception("Enzyme: multiple '%s'", enz->name.c_str() );
}


void  Cell:: computeRate( Solution &dS, Real t, Real E, const Solution &S_in )
{
	const Real      zeta_m = E * F_over_RT;
	const Solution &S_out  = solutions[ iS_out ]; //! here: constant
	
	//--------------------------------------------------------------------------
	//-- initialize dS
	//--------------------------------------------------------------------------
	dS.zero();

	//--------------------------------------------------------------------------
	//-- passive fluxes
	//--------------------------------------------------------------------------

	//(std::cerr << "l").flush();
	dS.addLeak(zeta_m, Volume, S_in, S_out);
	//std::cerr << "dS_leak=" << std::endl; dS.display(std::cerr);
	
	//--------------------------------------------------------------------------
	//-- enzymatic fluxes
	//--------------------------------------------------------------------------
	//(std::cerr << "e").flush();
	for( EnzymeDB::iterator i = enzymes.begin(); i != enzymes.end(); ++i )
	{
		EnzymePtr &enz = *i;
		(void) enz->addRate(dS,t,zeta_m,S_in,S_out,*this);
	}
	//std::cerr << "dS_leak+enz=" << std::endl; dS.display(std::cerr);

	//--------------------------------------------------------------------------
	//-- chemical reactions
	//--------------------------------------------------------------------------
	// (std::cerr << "u").flush();
	chemsys.updateRate( dS, t, S_in );
	
}

void  Cell:: computeDiff( array<Real> &dFdt, Real t, const array<Real> &F )
{
	assert( F.size() >= nvar );
	assert( dFdt.size() == F.size() );
	Solution &S  = solutions[ iS_tmp ];
	Solution &dS = solutions[ idS ];
	const size_t jEm = M+iEm;
	

	//--------------------------------------------------------------------------
	//-- Extract variables
	//--------------------------------------------------------------------------
	S.load( F );
	const Real E = F[ jEm ];
	
	//--------------------------------------------------------------------------
	//-- Compute the net increase with current potential and concentrations
	//--------------------------------------------------------------------------
	computeRate( dS, t, E, S );
	
	//--------------------------------------------------------------------------
	//-- Compute the net charge change: 
	//-- the law of mass action doesn't change it !
	//-- Volume is initially in microns^3 => L
	//--------------------------------------------------------------------------
	const Real Im  = 1.0e-15 * Faraday * Volume * dS.sumCharge(); 

	//--------------------------------------------------------------------------
	//-- store solution information
	//--------------------------------------------------------------------------
	dS.save(dFdt);
	
	//--------------------------------------------------------------------------
	//-- store extra information
	//--------------------------------------------------------------------------
	dFdt[ jEm ] = Im/Capa;
	
}

Real Cell:: DeltaQ() const
{
    return (solutions[ iS_in].sumCharge() - solutions[ iS_out ].sumCharge()) * ( 1.0e-15 * Faraday * Volume );
}

Real Cell:: DeltaEm() const
{
    return DeltaQ()/Capa;
}

Real Cell:: excessConc( const Real Erest) const
{
    const Real Q_real = Erest*Capa;
    const Real dQ     = Q_real - DeltaQ();
    return dQ/( 1.0e-15 * Faraday * Volume );
}

